A mechanohydraulic model supports a role for plasmodesmata in cotton fiber elongation

Abstract Plant cell growth depends on turgor pressure, the cell hydrodynamic pressure, which drives expansion of the extracellular matrix (the cell wall). Turgor pressure regulation depends on several physical, chemical, and biological factors, including vacuolar invertases, which modulate osmotic pressure of the cell, aquaporins, which determine the permeability of the plasma membrane to water, cell wall remodeling factors, which determine cell wall extensibility (inverse of effective viscosity), and plasmodesmata, which are membrane-lined channels that allow free movement of water and solutes between cytoplasms of neighboring cells, like gap junctions in animals. Plasmodesmata permeability varies during plant development and experimental studies have correlated changes in the permeability of plasmodesmal channels to turgor pressure variations. Here, we study the role of plasmodesmal permeability in cotton fiber growth, a type of cell that increases in length by at least three orders of magnitude in a few weeks. We incorporated plasmodesma-dependent movement of water and solutes into a classical model of plant cell expansion. We performed a sensitivity analysis to changes in values of model parameters and found that plasmodesmal permeability is among the most important factors for building up turgor pressure and expanding cotton fibers. Moreover, we found that nonmonotonic behaviors of turgor pressure that have been reported previously in cotton fibers cannot be recovered without accounting for dynamic changes of the parameters used in the model. Altogether, our results suggest an important role for plasmodesmal permeability in the regulation of turgor pressure.


Model parameters
Values for the model parameters listed in the main text in Table 1 are determined based on a literature review.In the absence in the literature of direct values for the plasmodesmal permeability and the solute source in cotton fibers, we estimated here these two parameter values, using other measurements from the literature.Parameter values used for these estimations together with the source references are summarised in Table S1.More precisely, for each parameter we identified a range of observed values, and we compute a corresponding reference value as the geometric mean of the range limits.We also verified that the typical values used in other examples of modeling approaches to plasmodesmal permeability, like in (1,2), lie within the range that we used.

A. Estimation of plasmodesmal permeability.
To estimate plasmodesmal permeability µ, we followed the methodology of (3) and first calculated the conductivity of a single plasmodesma using different laws depending on the assumptions that we made, see Fig. S1.
If a plasmodesma were just an empty cylindrical tube of radius R and length w, its permeability λ R to a fluid of viscosity η would be given by the Hagen-Poiseuille law, λ R = πR 4 8ηw .[1] However, the desmotubule obstructs the middle of the plasmodesma and can be considered as a cylinder of radius γR with γ ∈ [0, 1].Therefore, even if we assume that the space between the desmotubule and the plasma membrane (cytoplasmic sleeve) is free, the Hagen-Poiseuille law has to be modified to account for the occlusion of the plasmodesma channels by the desmotubule, see Figures S1-a for the transversal and S1-b for the longitudinal section.In this case, the equation for annular flow yields the permeability of a single plasmodesma, (1,4,5) We plotted the dependence of this permeability on the relative desmotubule radius γ in Figure S1-g.Notice that the permeability drastically drops with the presence of the desmotubule: it loses one order of magnitude with respect to the empty cylinder with a relative desmotubule radius of 0.5.For desmotubules that are not very thin, this permeability can be approximated quite well (see Figure S1-g) by the following polynomial expression in γ λ central_proxy = 2 3 λ R (1 + γ) (1 − γ) 3 . [3] If we consider now that the space between the desmotubule and the plasma membrane is obstructed so that the cylindrical channels of radius k in the ultrastructure of the plasmodesma take half of the annular volume left by the desmotubule (6)(7)(8), see Figure S1-c, then the permeability of a single plasmodesma becomes Notice that the channels have to fit in the cytoplasmic sleeve, therefore, they have a maximal radius k max = R(1 − γ)/2, given the plasmodesma and desmotubule radii.In Figure S1-g, we represented the permeability λ obstructed as a function of γ with the reference value for k based on measurements (see Table S1), unless this exceeds the maximal value k max for a given γ value, where we considered k = k max instead.Figure S1-g allows us to compare the predictions for the permeability of a single cylindrical plasmodesma given the above-presented models: central desmotubule with and without obstructed cytoplasmic sleeve.In particular, for the same radius R, length w, and relative desmotubule radius γ, we can get an upper estimation for the permeability by considering the empty cytoplasmic sleeve, as well as a lower estimation by considering an obstructed plasmodesma with channels.Therefore, by considering the whole biologically relevant parameter space (see Table S1) for the parameters appearing in these two models (central and obstructed), and taking the minimal and the maximal plasmodesmata permeability value, we estimate a range for the plasmodesmata permeability λ.
We note that other studies have explored more complex geometries, such as a desmotubule with an offset or varying radius (2), necks (9) or funnel-shaped plasmodesmata (10).Such geometries lead to changes in the permeability of a plasmodesma with respect to the above-presented geometries based on an ideal cylinder.In the following, we briefly present some of these more complex geometries as perturbations of the cylindric one, so that we can estimate the effect of these geometrical perturbations on the permeability of one plasmodesma.
On the one hand, in (2), they study for instance the effect of a possible offset of the desmotubule inside the plasmodesma.Let δ ∈ [0, 1] be the offset parameter, namely the distance between the axes of the desmotubule and of the plasmodesma is δ(1 − γ)R, as schematized in S1-c.Then the permeability of a single plasmodesma with a relatively wide desmotubule increases quadratically with the offset parameter δ In particular, the permeability increases as soon as there is an offset of the desmotubule, and it can increase by a maximal value of 150% with respect to the unperturbed concentric plasmodesma when the desmotubule is in contact with the plasma membrane (δ = 1).We traced the permeability λ of f set (δ = 1) in this extreme situation in figure S1-h.
On the other hand, a neck reduces permeability in the same proportions as an obstruction.To see this, we consider a plasmodesma with a neck region as presented in figure S1-e.In this context, the neck represents an annular obstruction of thickness ω(1 − γ)R with ω ∈ [0, 1] and length ζw with ζ ∈ [0, 1].This structure can be seen as a series of two annular pipes with different outer radii, and therefore, we approximate its permeability by We represented in figure S1-h the permeability corresponding to this model taking into account a neck region which obstructs 50% of the cytoplasmic sleeve thickness (ω = 0.5), on 10% of the length of the plasmodesma (ζ = 0.1).One can see that indeed, the presence of the neck region reduces the permeability of the plasmodesma, and its effect is of the same order of magnitude as of the obstruction in the obstructed model.
Furthermore, the widening on one end of funnel-shaped plasmodesmata can drastically raise their permeability.Considering a plasmodesma of radius R on the one end and of radius (1 + )R on the other end, with a concentric desmotubule of radius γR, as shown in figure S1-f, we obtain the following expression for its permeability We represented in figure S1-h this function of the relative desmotubule radius γ for the specific values of = 0.5 for a 50% and = 1 for a 100% raising of the plasmodesmata radius on the one end.We see that an order of magnitude in permeability can easily be gained from this type of widening.However, since we did not find in the literature evidence for the existence of these funnel-shaped plasmodesmata for cotton on the cell wall separating the seed and the fiber, we did not consider this model in the estimation of the biological range of plasmodesma permeability.
Finally, the value of the conductivity of a plasmodesma, λ, was then scaled to the cell by multiplying it by the surface area of the cell wall shared by the fiber and the contiguous epidermal cell, and ρ the surface density of the plasmodesmata.This yields the plasmodesmal permeability of the fiber, Fig. S1.The effect of the geometry on the permeability of plasmodesmata.Different models of a plasmodesma are presented as follows: a) transversal section of a plasmodesma with a concentric desmotubule; b) longitudinal section of a cylindrical plasmodesma with concentrical desmotubule; c) transversal section of an obstructed plasmodesma; d) transversal section of a plasmodesma with a desmotubule with an offset; e) longitudinal section when necks obstruct the plasmodesma; f) longitudinal section of funnel shaped plasmodesmata.g) Permeability λ as a function of the relative desmotubule radius γ.The effect of an offset desmotubule and of the obstruction of the cytoplasmic sleeve in the case of a cylindrical plasmodesma are shown.h) The effect of necks and of a funnel shape are compared to the effect of the preceding perturbations.For g) and h) unless explicitly mentioned, all parameter values take the reference values from Table S1.The reference value for the relative desmotubule radius (γ = 0.82) is indicated by a vertical black line on the graphics.Calculated with equation ( 8), geometric mean of the maximum and minimum permeability values enabled by the parameter's range.
Recall, that in order to estimate the range for the total plasmodesmal permeability µ in Table S1, we considered all the parameter space defined by the ranges for the other parameters in this table, considered the models of non-obstructed (equation 2) and obstructed (equation 4) cylindrical plasmodesma with concentrical desmotubule, and took the minimal and maximal value for µ.We note that if we would have taken just the reference value of the parameters to estimate the range of µ, then we would have obtained the narrow interval 21.2 × 10 −24 − 87.9 × 10 −24 m 3 P a −1 s −1 , and a reference value of 43.1 × 10 −24 m 3 P a −1 s −1 which is very close to the value obtained in Table S1, computed using the whole parameter space.

B. Estimation of the source of solute.
Since we assumed that π fiber is proportional to the concentration of solutes, we estimated the typical value of the source parameter, α, as π fiber κ, where π fiber was extracted from (17) and κ is the average value of growth rate that we estimated from the increase of fiber length with time.We extracted experimental data of the length of cotton fiber cells reported in (18)(19)(20)(21).We fitted lengths to the logistic growth equation: where M is the maximum volume and t 0 is the midpoint of the curve.We then used the average value of κ.

Analytical study of the behaviors of the system
Here we aim to mathematically address whether the model with constant parameters can predict a peak in (turgor or osmotic) pressure.To do so, we non-dimensionalized the governing equations and we studied the sign of the second derivative of the turgor and osmotic pressure with respect to time.We considered all possible cases, according to dimensionless parameters and we analytically derived conditions for the existence of a peak in pressure.This section contains the derivation of the equations to identify the different possible behaviors of the system.The main question is to know whether the model can have non-monotonic behaviors for the turgor and osmotic pressure like reported in experimental data.
Throughout this analysis we use the following dimensionless form of the equations of the model The dimensionless variables are defined from turgor pressure, osmotic pressure, cell volume, as And the dimensionless parameters are defined by Table S2.List of the range of the dimensionless parameters with parameters values of table 1.

Symbol Expression Range
To make it easier for the reader, we only provide here the analysis of the behavior of p.We note that we obtained the same results for π.
A. Case without growth.When turgor pressure in the fiber P fiber is lower than the yield threshold Y , we get p = vπ v + 1 .
Replacing in equation ( 10) we obtain for the volume: The volume remains constant, which closes this case.

B. Case with growth.
Throughout this part, we assume that turgor pressure in the fiber P fiber is higher than the yield threshold Y (in terms of the dimensionless parameters p > y 0 ).We consider the equations ( 10) and ( 11) and the dimensionless turgor pressure p is It is easy to check that the first derivative of v is always positive if p > y 0 .From Eqs. ( 20) and ( 10), we obtain dv dτ = vφ 0 (p − y 0 ) > 0. [21] The first derivative of the dimensionless turgor pressure, p, is Because Eq. ( 11) contains the positive part of p we distinguish the cases when p is strictly negative (section B.1) and when p is positive or null (section B.2).
First of all, we write below the second derivative of the dimensionless volume, v, the dimensionless osmotic pressure, π, and the dimensionless turgor pressure, p, that we will use in the following sections.Equation ( 10)) yields Finally, Equation (22) results in In this analysis, we seek whether p reaches a maximum.We thus assume the existence of a dimensionless time τ 0 ≥ 0 such that the first derivative of the dimensionless turgor pressure p is equal to zero dp dτ (τ 0 ) = 0.
This can be rewritten using Eq. ( 22) as Then, using the Eq. ( 25), we find a necessary condition to have the second derivative of the dimensionless turgor pressure p negative Hereafter we consider p and other quantities to be evaluated at time τ 0 , dropping (τ 0 ) for brievity.By developing the equation (26) using the differential equations ( 10) and ( 11) and the expression of p (20) we obtain We obtain therefore a quadratic polynomial of p.This polynomial can have either zero, one or two solutions corresponding to the value of the dimensionless turgor pressure p at dimensionless time τ 0 when the first derivative of p equals zero.We study then independently the different cases depending on the sign of p. B.1.Case P fiber < P seed (p < 0).We seek to demonstrate in this section that the dimensionless turgor pressure p → p(τ ) defined over R + can have at least one local maximum, implying that the function p can be non-monotonic with a local maximum.Because p < 0 in this section, the dimensionless yield threshold y 0 must be strictly negative, y 0 < 0. We reformulate Eq. (27) using p < 0. At τ = τ 0 we have The discriminant of this quadratic polynomial is Let us study the sign of the discriminant ∆.The discriminant is positive if and only if Since y 0 < 0, all terms of this quadratic polynomial of v (in Eq. 29) are strictly positive.Therefore ∆ > 0 and Eq.(28) has two solutions p + and p − such that .
So far, we have analyzed whether the first derivative of p(τ ) vanishes.We found that the only solution satisfying p (τ 0 ) < 0 and p (τ 0 ) > y 0 is p + .We next study the sign of the second derivative of the dimensionless turgor pressure p.As we assume that p < 0, Eq. (24) becomes For τ = τ 0 (implying that dp/dτ = 0), this yields Meanwhile, Eq. ( 23) becomes Then, replacing the expression of the two terms above in equation Eq. ( 22) we obtain We use Eq. ( 26), and we obtain We recognize the equation Eq. ( 10) in the parenthesis, therefore this expression can be simplified to Since we assume p < 0 in this section, we finally have the following equivalence: Under the assumption p > y 0 , the first derivative of the dimensionless volume v is strictly positive.Thus we can only have a non-monotonic behavior with a local maximum (or a monotonic behavior if the first derivative does not vanish) and this is possible only if α 0 < φ 0 y 0 (φ 0 y 0 − 1).
We notice that if α 0 ≥ φ 0 y 0 (φ 0 y 0 − 1) then the first derivative of p is strictly positive when p is in the interval [p − , 0].Because p − < y 0 , we must have p > p − in this case of growth.As a consequence, turgor pressure in the fiber increases until it reaches the value of the turgor pressure in the seed cells, thus the case p < 0 is no longer true and we get p ≥ 0, which is the case investigated in next section B.2.
Repeating the analysis yields the same results for the dimensionless osmotic pressure π -we do not present this analysis here.B.2. Case P fiber ≥ P seed (p ≥ 0).We seek to demonstrate in this section that the dimensionless turgor pressure p → p(τ ) defined over R + can have at most one local minimum but no local maximum.Because we assume that p ≥ 0, the dimensionless yield threshold y 0 can be either positive or negative.
We rewrite Eq. ( 27) using p ≥ 0. At τ = τ 0 we have [33] The discriminant of this last quadratic polynomial is There are solutions to Eq. 33 only if ∆ ≥ 0, which is equivalent to V. Hernández-Hernández, O. C. Marchand, A. Kiss, and A. Boudaoud This is a quadratic polynomial for v, with discriminant ∆ .If ∆ ≤ 0, then condition (34) is true because the coefficient of v 2 is strictly positive.If ∆ > 0, then the polynomial in v has two roots v 1 and v 2 (v 1 < v 2 ).If both v 1 and v 2 are negative, then, because we consider v > 0, condition (34) is true.If only one solution is positive (it would be v 2 ), then we must have v ≥ v 2 for Delta ≥ 0. If both solutions are positive, then we must have v ≤ v 1 or v ≥ v 2 for Delta ≥ 0. Now, we have all the conditions to have ∆ ≥ 0. In the following, we assume the inequality to be strict as the case ∆ = 0 can be obtained by continuation.The roots of the polynomial in Eq. 33 are p − and p + , such that We next examine whether these roots meet the conditions of this section, i.e. p (τ 0 ) ≥ 0 and p (τ 0 ) > y 0 .We start with The outcome depends on the sign of g(v).If negative, then the solution p − does not meet the condition p − ≥ 0. For the solution p + we obtain is positive, then the solution p + satisfies the condition, and for the other solution we have p − (τ 0 ) ≥ 0 ⇐⇒ α 0 ≤ φ 0 y 0 (φ 0 y 0 − 1) .Now, we consider the condition p (τ 0 ) > y 0 .We have If f (v) > 0, then the solution p − does not satisfy the relation p (τ 0 ) > y 0 , and for the solution p + we obtain If f (v) < 0, then the solution p + satisfies the relation p (τ 0 ) > y 0 , and for the solution p − we have after calculation We remark that Thus, we can not have simultaneously for the same dimensionless volume v the condition f (v) < 0 and g (v) > 0, thus the solution p − does not meet simultaneously the conditions p − (τ 0 ) > y 0 and p − (τ 0 ) > 0. The only solution for the case ∆ > 0 to have the first derivative equal to zero is therefore p + , under some conditions on v according to the value of the parameters.
In what follows we study the sign of the second derivative of p. Taking into account that p ≥ 0, the second derivative of the dimensionless osmotic pressure π (Eq.24) can be rewritten as The second derivative of the dimensionless volume v (Eq.23) becomes At τ = τ 0 , these equations become We can calculate the second derivative of p (Eq. 22), and compare its value to zero The case p = 0 does not meet this condition.Because we assumed that the fiber is growing, dv/dτ (τ 0 ) > 0, thus the inequality above implies that Because we assumed that p is positive, the quantity φ 0 y 0 − 1 must be positive, we therefore have φ 0 y 0 > 1 which also implies that y 0 > 0 because φ 0 > 0. The solution p + must verify Eq. (37) for the second derivative of p to be negative, which is equivalent to If q (v) < 0, the p + (τ 0 ) does not verify inequality (37).If q (v) > 0 then we obtain We assumed that φ 0 y 0 > 1 thus y 0 > 0. We thus find that Since y 0 > 0 and v > 0 we have the following implication But we assumed that φ 0 y 0 > 1 for (37), and Thus we obtain a contradiction, and the case q (v) > 0 is not compatible with Eq. (37).Thus the solution p + does not verify the condition (37), and the second derivative of p can not be negative.
A similar analysis shows that it is possible to have the second derivative of p positive with the solution p + .

B.3. Relations with the initial conditions.
Here we consider whether the initial conditions enable an extremum of pressure.In the preceding section, we investigated the sign of the second derivative of p.We obtained the value of p when the first derivative is equal to zero (Eq.35).We can also calculate the isocline of p in the phase space (p, v) and show that this isocline is strictly monotonic, increasing for p ≥ 0 and decreasing for p < 0. This leads to a necessary condition for the first derivative to vanish at a certain time τ 0 .Indeed, for p to have a local maximum the initial condition p(0) must be lower than the maximum reached, and for the local minimum, the initial condition must be higher than the minimum reached (i.e. when the first derivative equals to zero).We note that we also obtained the same results for π.
In the case p < 0, in order to have the second derivative of the turgor pressure negative, it is necessary to have p (0) < p max such that with ∆ defined as In the case p ≥ 0, in order to have the second derivative of the turgor pressure positive, it is necessary to have p (0) > p min such that with ∆ defined as For the same reasons, in order to have the second derivative of the osmotic pressure negative, for the case p < 0, it is necessary to have π (0) < π max such that with ∆ defined as Finally, in order to have the second derivative of the osmotic pressure positive, for the case p ≥ 0, it is necessary to have π (0) > π min such that with ∆ defined as ∆ = 4α 0 (φ 0 v + 1) (v (φ 0 + 1) + 1) + (vφ 0 (y 0 + 1) + 1) 2 .
Moreover, because the osmotic pressure and the turgor pressure must be positive to have a physical meaning, we must have π ≥ −1, which implies p ≥ −v v+1 if p > y 0 or p ≥ v(φ 0 y 0 −1) v(1+φ 0 )+1 if p ≤ y 0 and P fiber > 0 implies p ≥ − P seed π seed .This naturally adds conditions on the parameters.For instance, for the minimum possible value of π, i.e. π = −1, to have the system has the osmotic pressure and the turgor pressure positive, we must have: or in other terms: We can deduce such relations between the parameters by considering the limit of p and π at +∞.These relations have to be verified if we want to have the turgor pressure and osmotic pressure positive in numerical simulations.
At this point, we have analysed all possible behaviors in terms of temporal variations of the three dimensionless variables v, π, and p.These behaviors and the corresponding conditions are summarised in the following section.

Summary of the analytical study and numerical illustrations
The different possible behaviors for the pressures and the volume with the model are detailed in the following figures and tables.
Table S3.Behaviors of the model.the quantities p max , p min , π max , π min are defined previously in equations (38,39,40,41) and depend on the initial dimensionless volume v (0) and on the three dimensionless parameters α 0 , φ 0 and y 0 .The variables v, p and π represent respectively the dimensionless volume, turgor pressure and osmotic pressure.Table S4.Values used for the curves in Figure S2.All the cases from Table S3 are illustrated.In all cases, we fixed the following parameters: L r = 0.1 M P a −1 .h−1 , µ = 0.001 mm 3 M P a −1 .h−1 , π seed = 1 M P a, P seed = 0.18 M P a, Y = 0.06 M P a, and φ = 1.2 M P a −1 .h−1 .Therefore φ 0 y 0 (φ 0 y 0 − 1) = 3.51 and 1 − y 0 φ 0 /φ 0 = 0.20; we state in this Table whether these quantities are higher or lower than α 0 .In the fist column, P fiber or nd π fiber indicates that, respectively, turgor or osmotic pressures are shown.α is given in units of M P a.h −1 , while the initial condition P (0) is in M P a.   S3.The values of the parameters used are given in Table S4.

Condition for a peak of turgor pressure and biologically-relevant parameters
We demonstrated that, in order to have a peak of turgor pressure, we need to have α 0 < φ 0 y 0 (φ 0 y 0 − 1) and y 0 < 0. In figure S3 we plot this condition in the dimensionless parameter space using parameter values reported in table S2.We note that the values of φ 0 y 0 (φ 0 y 0 − 1) are very low compared to the biophysical value reported for α 0 .The necessary conditions α 0 < φ 0 y 0 (φ 0 y 0 − 1) and y 0 < 0 can be satisfied if the value of α 0 is lower than approximately 6.0 × 10 −4 , which represents 1.7% of the biophysical range reported of α 0 .S1, and the continuous purple line corresponds to this value of α = α 0 min .

Details of the sensitivity analysis
As part of the sensitivity analysis of the model near reference parameters, we show in Figures S4-S6 the ratio between the final value of the observable variable and its value for the reference parameters.

Dynamic source of solute
In Figure S7 we plot the behavior of the volume and pressures with a dynamic source of solute α corresponding to the temporal pattern shown in Panel a.

Fig. S2 .
Fig. S2.Behaviors of the osmotic and turgor pressure according to the cases reported in TableS3.The values of the parameters used are given in TableS4.

Fig. S3 .
Fig. S3.Condition α 0 ≥ φ 0 y 0 (φ 0 y 0 − 1) in the space of dimensionless parameters φ 0 and y 0 .The hatched area corresponds to the region of the parameter space where the condition is satisfied, for α 0 min = 1.0 × 10 −5 as computed from TableS1, and the continuous purple line corresponds to this value of α = α 0 min .

Fig. S4 .Fig. S5 .Fig. S6 .
Fig. S4.Sensitivity of the volume to the seven parameters of the model.Normalised value plotted as function of the relative change of the parameter.

Fig. S7 .
Fig. S7.Pressures and volume with a dynamic source of solute α. a) Source of solute α as a function of time.b) Normalised final volume as function of the start time t 0 and of the duration t dist of the period during which the source is enhanced (the normalised volume corresponds to the ratio of the final volume at 3000 h to the final volume at the same time taking t dist = 50 h and t 0 = 200 h)).c) Continuous lines correspond to the behavior of the model with dynamic source.Dashed lines correspond to the case when all parameters are constant and equal to their reference values.